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Abstract. This article focuses on automatically generating polynomial 
equations that are inductive loop invariants of computer programs. We 
propose a new algorithm for this task, which is based on polynomial 
interpolation. Though the proposed algorithm is not complete, it is ef- 
ficient and can be applied to a broader range of problems compared to 
existing methods targeting similar problems. The efficiency of our ap- 
proach is testified by experiments on a large collection of programs. The 
current implementation of our method is based on dense interpolation, 
for which a total degree bound is needed. On the theoretical front, we 
study the degree and dimension of the invariant ideal of loops which have 
no branches and where the assignments define a P-solvable recurrence. 
In addition, we obtain sufficient conditions for non-trivial polynomial 
equation invariants to exist (resp. not to exist). 



1 Introduction 

Many researchers have been using computer algebra to compute polynomial loop 
invariants, see for instance [21,22,19,17,14,15,9,3,1,20,5,12]. In this article, 
we propose an alternative method, based on interpolating polynomials at finitely 
many points on the reachable set of the loop under study. This interpolation 
process 3 yields "candidate loop invariants" which are checked by a new criterion 
based on polynomial ideal membership testing. 

Our paper proposes the following original results. On the theoretical front, for 
P-solvable loops with no branches, we supply a sharp degree bound (Theorem 1) 
for the invariant ideal, as well as dimension analysis (Theorem 3), of the invariant 
ideal. We establish a new criterion (Corollary 1) based on polynomial system 
solving for checking whether or not a given conjunction of polynomial equations 
is indeed a loop invariant. Meanwhile, Corollary 2 states a sufficient condition 
for the invariant ideal of a loop to be trivial. 

On the algorithmic front, we propose a modular method (Algorithm 2) for 
generating polynomial loop invariants. Thanks to polynomial interpolation, most 
of our calculations reduce to linear algebra. As a consequence, the proposed 
method works in time n d , where n is the number of loop variables and d is 
the total degree of polynomials to interpolate. 



3 Note that polynomial interpolation is different from the interpolation in [13], which 
is called Craig interpolation in first order logic. 



Our method is probabilistic and may not compute the whole invariant ideal. 
However, the implementation (in Maple) of our method computes all the in- 
variants given in each example proposed by Enric Rodriguez Carbonell on his 
page 4 . Moreover, the degree and dimension estimates can help certifying that 
whole invariant ideal has been obtained. For instance in co-dimension one, the 
invariant ideal is necessarily principal. 

Our method needs not to solve the recurrence relations associated with the 
loops and thus does not need to manipulate the algebraic numbers arisen as 
eigenvalues of these recurrence relations. Therefore, all polynomials and matrices 
involved in our method have their coefficients in the base field. Our method 
applies to all loops which can be modeled as algebraic transition systems [21] 
and can be generalized to handle loops which can be modeled as semi-algebraic 
transition systems [3] . It can be applied to compute all kinds of invariants (see 
the notions of different loop invariants presented in Section 2), which is not the 
case for the methods based on "recurrence solving" [20, 12]. In particular, the 
methods in [20] and [12] apply only to compute absolute inductive invariants, 
that is, the loop guard and branch conditions are ignored (thus the loop goes 
to the branches randomly). This means that those methods can not find loop 
invariants which are not absolute inductive invariants. 

Our implementation is tested against the implementation of 3 other meth- 
ods [20, 5, 12] which were kindly made available to us by their authors. The exper- 
imental results are shown in Section 4. While the performance of our method is 
comparable to that of [5] on the tested (somehow simple) examples, our method 
has less restrictive specifications than the methods of [20, 12] which target only 
on absolute polynomial loop invariants for P-solvable loops/solvable mappings. 
In addition, the method in [20] applies only when all assignments are invertiblc 
and the eigenvalues of coefficient matrix for the linear part are positive ratio- 
nals, and it is claimed to be complete, that is, to compute all absolute polynomial 
loop invariants, when it is applicable. The method in [12] can be applied all P- 
solvable loops in theory and is complete for loops without branches. However, in 
the implementation, the assignment in the loop are required to be not "coupled" 
together. 

Let us conclude this introduction with a brief review of other works on loop 
invariant computation. In [10], linear equations as invariants of a linear program 
at each location is considered, by tracking the reachable states with a method 
based on linear algebra. In [15], the method in [10] is improved and generalized 
to generate polynomial equations as invariants; it is also shown there that check- 
ing whether or not a linear equation is invariant is undecidable in general. In 
[14], for polynomial programs, the Authors discuss methods based on abstract 
interpretation, on checking whether or not a given polynomial equation is in- 
variant, as well as generating all polynomial invariants to a given total degree. 
In [17,5], a different abstract interpretation technique is developed, which uses 
polynomial ideal operations (e.g. intersection, quotient) as widen operators. In 
[9] and [3], quantifier elimination techniques are used to infer invariants from a 

4 http : / /www . lsi . upc . edu/~ erodri/webpage/polynomial_invar iant s/list . html 



given template; these methods requires expertise on supplying meaningful tem- 
plates, while the complexity of quantifier elimination also restrict their practical 
efficiency. In [21], Grobner basis together with linear constraint solving is used 
to infer polynomial equations as invariants. 

2 Preliminaries 

Let Q denote the rational numbers and Q the algebraic closure of Q. Let Q* 
(resp. Q ) denote the non zero elements in Q (resp. Q) . 

2.1 Notions on loop and loop invariants 

We will use the following simple loop (in MAPLE-like syntax) to introduce some 
notions related to loop and loop invariant that we are going to use. 



x := a; 




V ■= b; 




while x 


< 10 do 


x := 


= x + y 5 ; 


y ■= 


= i/ + i; 


end do; 





A loop variable of a loop is a variable that is either updated in the loop; 
or used to initialize/update the values of other loop variables, e.g. x,y,a,b are 
loop variables. Without loss of generality, we assume that all variables take only 
rational number values, i.e. from Q. By initial values of a loop, we mean all 
possible tuples of the loop variables before executing the loop; the set of the 
initial values of the above loop is 

{(x, y,a 7 b) \x = a,y = b, (a, b) e Q 2 }. 

Given an initial value v, the trajectory of the loop starting at v, is the sequence of 
all tuples of of loop variable values at each entry of the loop during the execution, 
with the loop variable being initialized by v; the trajectory of the above loop 
starting at (x, y, a, b) = (1, 0, 1, 0) is 

(1,0,1,0), (1,1,1,0), (2,2,1,0), (34,3, 1,0). 

The collection of value tuples of all trajectories is called reachable set of the loop. 
Note that, in general, it is hard to describe a reachable set of a loop precisely. 
A loop invariant (or plain loop invariant) of a loop is a condition on the loop 
variables satisfied by all the values in the reachable set of the loop. 

By inductive reachable set of a loop, we mean the reachable set of the loop 
while ignoring the guard condition, while by absolute reachable set of a loop, 



we mean the reachable set of the loop while ignoring the guard conditions, the 
branch conditions and viewing branches to be selected randomly. Then, by an 
inductive (loop) invariant (resp. absolute (loop) invariant) of a loop is a condition 
on the loop variables satisfied by all the tuple values in the inductive (resp. 
absolute) reachable set of the loop. 

It is easy to deduce that an absolute invariant is always an inductive invari- 
ant, and an inductive invariant is always a loop invariant. In principle, absolute 
inductive invariants are easier to study and compute than the inductive invari- 
ants and plain invariants. However, the absolute invariants can be trivial, which 
is not of practical interest to program analysis. See the following example [21] 
for the case of a trivial absolute invariant while inductive loop invariants in not 
trivial. 



vi ■= 0; 




2/2 ■= 0; 




2/3 := ari; 




while t/2 7^ 


do 


if 2/2 + 1 


= x 2 


then 




2/1 = 


= 2/1 + 1; 


2/2 ■■ 


= 0; 


2/3 : 


= 2/3-1; 


else 




2/2 : 


= 272 + 1; 


2/3 : 


= 2/3-1; 


end if 




end do 





Indeed, the condition 1/1X2 + 2/2 + 2/3 = x i is an inductive invariant of the above 
loop. Note there are also loop invariants which are not inductive invariants, e.g. 
x — 1 = is an invariant but not an inductive of the following loop. 



x:=l; 




while x 


^ 1 do 


x := 


= x + 1; 


end do 





On the other hand, the inductive invariants are less likely to be trivial and 
easier to handle than the loop invariants. In this article, we are interested in 
the inductive invariants that are given by polynomial equations and that we call 
polynomial equation invariants, or simply polynomial invariants when there is 
no possible confusion. It is not hard to deduce that all polynomials that are 
inductive invariants (or loop invariants, or absolute invariants) of a loop form 



an ideal (which is indeed the ideal of the points in the inductive reachable set), 
one can also refer [18] for an alternative proof. We call the ideal of polynomials 
which are inductive invariants of a given loop the invariant ideal of the loop. 

In this paper, we consider loops of the following shape. 



while Co do 




if d 




then 




X 


= A 1 (X); 


elif C 2 




then 




X 


= A 2 (X); 


elif C m 




then 




X 


= A m (X); 


end if 




end do 





where 

1. X = x\, X2, ■ ■ ■ ,x s is a list of s scalar loop variables, taking values from Q; 

2. the initial values of the loop are constrained by polynomial equations and 
polynomial inequations; 

3. the Cj's are pairwise exclusive algebraic conditions (polynomial equations 
and polynomial inequations) on X; 

4. the Ai's are polynomial functions of X with coefficients from Q. 

In our loop model, when the loop body contains assignments only (thus no 
branches), the assignment indeed induces a recurrence relation among the loop 
variables, which are viewed as recurrence variables. In this case, we shall simply 
refer to the loop as this recurrence relation and the initial values of the loop. 
Here we will show briefly how the loop invariant can be computed by explicitly 
solving the recurrence relation. Later, in our theoretical analysis, presented in 
Sections 3, we will focus on the study degree and dimension of the invariant 
ideal of such kind of loops, where the induced recurrence relation is so-called 
P-solvable recurrence. 

Example 1 Consider the loop computing the sequence of the Fibonacci numbers: 



x := 0; 

while true do 

(x,y) := (y,x + y); 
end while 



Viewing (x,y) as two recurrences variables, the loop is actually computing the 
two recurrence sequences of values of x and y defined by the following recurrence 
relation and initial condition: 



x(n + 1) = y(n), y(n + 1) = x(n) + y(n), with x(0) = 0, y(0) = 1. 
We can write down the closed form for x(n) and y(n) as follows: 

( y/5 + l ^ I -V5 + 1 ^ 



V5 Vb 

5+1 i^r - 

2 ^5 2 v/5 ' 

Let a,u,v be 3 variables. Replace ( (Vesp. ( ~ V ^ +1 )"J fry u (Vesp. % wj; 

replace \/5 &?/ a. Taking the dependencies u 2 v 2 = 1, a 2 = 5 on f/ie new variables 
into account, the invariant ideal of the loop is 

, au av a + lu —a+lv , _ 9 9 

which turns out to be (1 — y 4 + 2.xy 3 + .x 2 y 2 — 2x 3 y — x 4 ). 
2.2 Poly-geometric summation 

As we discussed in the previous subsection, the study of loops without branches 
can be reduced to the study of recurrence sequences. In this subsection, wc recall 
several well-known notions together with related results adapted to our needs. 

Those notions and results can usually be stated in a more general context, 
e.g. the notion of multiplicative relation can be defined among elements of an 
arbitrary Abelian group, whereas we define it for a multiplicative group of alge- 
braic numbers. 

Definition 1 Let a\, . . . ,ctk be k elements of Q \ {1}. Let n be a variable tak- 
ing non-negative integer values. We regard n, a™, . . . , a£ as independent variables 
and we call a™, . . . , a£ n-exponential variables. Any polynomial o/Q[n, a™, . . . , a£" 
is called a poly-geometrical expression in n over Q w.r.t. ai, . . . , otk- 

Let f, g be two poly-geometrical expressions n over Q w.r.t. ai, . . . , au- Given 
a non-negative integer number i, we denote by f\ n =i the evaluation of f at i, 
which is obtained by substituting all occurrences of n by i in f. We say that f 
and g are equal whenever f\ n =i = g\ n =i holds for all non-negative integer i. 

We say that f{n) is in canonical form if there exist 



(i) finitely € Q , and 

(ii) finitely many pairwise different couples ei), . . . , (/3 m , e m ) all in (Q \ 

{1}) x Z>o, and 
(m) a polynomial Co(n) € Q[n], 

smc/i iftai eoc/i . . . , /3 m is a product of some of the a±, . . . , and smc/i i/iat ifte 
poly-geometrical expressions f(n) and^2^ =1 Ci /?" n ei + Co(n) are equal. When 
this holds, the polynomial Co(n) is called the exponential-free part of f(n). 

Remark 1 A^ote £/ia£ sometime when referring to poly-geometrical expressions, 
for simplicity, we allow n-exponential terms with base or 1, that is, terms with 
0™ or 1™ as factors. Such terms will always be evaluated to or 1 respectively. 

Proving the following result is routine. 

Lemma 1 With the notations of Definition 1. Let f a poly- geometrical expres- 
sion in n over <Q> w.r.t. on, . . . , ct!fc. There exists a unique poly- geometrical ex- 
pression c in n over Q w.r.t. ct\, ...,ah such that c is in canonical form and 
such that f and c are equal. We call c the canonical form of f. 

Example 2 The closed form f := ^ n+1 J " of X)"=o * 3 * s a poly- geometrical 
expression in n over Q without n-exponential variables. The expression g := 
n 2 2(" +1 ) — n 2™ 3^ is a poly- geometrical in n over Q w.r.t. 2, 3. Some evaluations 
are: /| (n=0 ) = 0,/|„=i = l,g\ n =o = 0,g\n=2 = 8. 

Notation 1 Let x be an arithmetic expression and let k € N. Following [6], we 
call k-th falling factorial of x and denote by x— the product 

x(x-l) ■■■(x-k + 1). 

We define x- := 0. For i = 1, . . . , k, we denote by {^} the number of ways to 
partition k into i non-zero summands, that is, the Stirling number of the second 
kind also denoted by S{n,k). We define {q} := 0. Finally, we shall make use of 
the convention 0° = 1. 

Example 3 The expression n 2 2^ n+1 ' — n 2™ 3^™/ 2 ' is clearly poly-geometrical in 
n over Q. Consider now a fixed non-negative integer k. The sum Y^i=i * fe has 
n — 1 terms while its closed form [6] below 

^ [if i + 1 

has a fixed number of terms and thus is poly-geometrical in n over Q. 

The following result is proved in [6]. 
Lemma 2 Let x be an algebraic expression and let k £N. Then we have 




Notation 2 Let r G Q and let k be a non-negative integer. We denote by 
H(r,k,n) the following symbolic summation 



H(r,k,n) :=J2 



r i-. 



i=0 



One can easily check that H(r,0,n) = r r _l holds for r ^ 1. Moreover, we 
have the following result. 

Lemma 3 Assume r^O. Then, we have 

(r-l)H{r,n,k) = {n - l) k r n - r k H(r, k - 1, n - 1). (1) 

Moreover, we have 

(i) if r = 1, then H{r,n,k) equals to ^"T/f ? which is a polynomial in n over Q 
of degree k + 1 . 

(ii) if r ^ 1, then H{r, n, k) has a closed form like r n f{n) + c, where f(n) is a 
polynomial in n over Q of degree k and c is a constant in Q. 

Proof: We can verify Relation (1) by expanding H(r,n,k) and H(r,k — 
l,n— 1). Now let us show the rest of the conclusion. First, assume r = 1. With 
Relation (1), we have 

kH(r,k-l,n-l) = (n-l) k . 

Therefore, we deduce 

n k±i 

F(r, »,*) = —. 

One can easily check that 2|== is a polynomial in n over Q and deg(s, n) = fc + 1. 

From now on assume r ^ 1. We proceed by induction on k. When k — 0, we 
have H (r, 0, n) = r r Zl ■ We rewrite T r Z\ as 

„ 1 1 



which is such a closed form. Assume there exists a closed form r" fk-i{n) +Ck-i 
for H(r, k — l,n), where fk-i{n) is a polynomial in n over Q of degree k — 1. 
Define 

_ (n - l)^r" - r fc (r"- 1 / fc _!(n - 1) + Cfc_i) 
r - 1 

It is easy to verify that s is a closed form of H (r, n, fc). We rewrite s as 

„ (n - l)-~ kfk-i{n - 1) _ rfcc fc _i 
r — 1 r — 1 

and one can check the later form satisfies the requirements of (ii) in the conclu- 
sion. This completes the proof. □ 



Lemma 4 Let k E N and let X be a non zero algebraic number over Q. Consider 
the symbolic summation 

n 
i=l 

1. if A = 1, fften f/iere exisis a closed form s{n) for S, where s is a polynomial 
in n over Q of degree k + 1 . 

2. if A ^ 1, £/ien t/iere errasfc a closed form A" s(n) + c /or S, w/iere s is a 
polynomial in n over Q of degree k and c e Q is a constant. 

Proof: By Lemma 2, we deduce 

ELi * fe ^ EIU fe} =1 {J}-') -v 
= Ej=i ({}} EI l =1 *a') 
= E" =1 ({, fe }^(A,j») 

Then, the conclusion follows from Lemma 3. □ 

The following definition of multiplicative relation specializes the general def- 
inition of multiplicative relation to non-zero algebraic numbers. 

Definition 2 (Multiplicative relation) Let k be a positive integer. Let A := 
(ai,...,afc) be a sequence of k non-zero algebraic numbers over Q and e := 
(ei, . . . , e/c) be a sequence ofk integers. We say that e is a multiplicative relation 
on A if n»=i a T = 1 holds. Such a multiplicative relation is said non-trivial if 
there exists i G {1, ...,n} such that e, ^ holds. If there exists a non-trivial 
multiplicative relation on A, then we say that A is multiplicativcly dependent; 
otherwise, we say that A is multiplicatively independent. 

All multiplicative relations of A form a lattice, called the multiplicative re- 
lation lattice on A, which can effectively be computed, for instance with the 
algorithm proposed by G. Ge in his PhD thesis [7]. 

For simplicity, we need the following generalized notion of multiplicative rela- 
tion ideal, which is defined for a sequence of algebraic numbers that may contain 
and repeat elements. 

Definition 3 Let A := (oti, . . . , a^) be a sequence of k algebraic numbers over 
Q. Assume w.l.o.g. that there exists an index i, with 1 < t < k, such that 
ati,...,cti are non-zero and a^+i, . . . , afe are all zero. We associate each on with 
a variable y ir where J/i,...,J/fc are different from each other. We call the mul- 
tiplicative relation ideal of A associated with variables y\, . . . ,yu, the binomial 
ideal of Q[yi, 2/2, • • • , Vk] generated by 

{ n v7- n vt vi 1 {vi,:-,v t )€Z} 

je{i,...J},vj>a ie{i,...,e},vi<o 



and {ye+i, ■ ■ ■ , Dk}, denoted by MRl(A; y\, . . . , yk), where Z is the multiplicative 
relation lattice on (a\, . . . ,otg). When no confusion is possible, we shall omit 
writ-ting down the associated variables yi,-.-,yk- 

Lemma 5 Let cti, . . . ,oik be k multiplicatively independent elements of Q and 
let n be a non-negative integer variable. Let f{n) be a poly-geometrical expression 
in n w.r.t. cti, . . . , ct^. Assume that f\( n =i) — holds for all i g N. Then, f is 
the zero polynomial of Q[n, a", . . . , aj}]. 

The following definition will be convenient in later statements. 

Definition 4 (Weakly multiplicative independence) Let A := (cti, . . . , a.k) 
be a sequence of k non-zero algebraic numbers over Q and let f3 € Q. We say (3 
is weakly multiplicatively independent w.r.t. A, if there exist no non-negative 
integers a, e 2 , . . . , such that (3 = n»=i a V holds. Furthermore, we say that 
A is weakly multiplicatively independent if 

(i) eti^l holds, and 

(ii) on is weakly multiplicatively independent w.r.t. 
{ai, . . . , (Xi-i, I}, for all i = 2, . . . , s. 

It is not hard to prove the following lemma on the shape of closed form 
solutions of single- variable linear recurrences involving poly-geometrical expres- 
sions. For the proof, we need the following lemma, which is easy to check, see 
for instance [16]. 

Lemma 6 Let n a variable holding non-negative integer values. Let a and b be 
two sequences in Q indexed by n. Consider the following recurrence equation of 
variable x: 

x(n) = a(n — 1) x(n — 1) + b(n — 1). 

Then we have 



x(n) = J] a(i) x(0) + 



i=0 



o rr s =o °( s ) 



Lemma 7 Let cti, . . . , a.^ be k elements in Q \ {1}. Let A e Q . Let h(n) be a 
poly-geometrical expression in n over Q w.r.t. ot\, . . . , a^. Consider the following 
single-variable recurrence relation R: 

x(n + 1) = Xx(n) + h(n). 

Then, there exists a poly-geometrical expression s(n) in n over Q w.r.t. cti, ■ ■ ■ , dk 
such that we have 

deg(s(n), a") < deg(/i(n), a") and deg(s(n), n) < deg(/i(n), n) + 1, 

and such that 



— if A = 1 holds, then s(n) solves R, 

— if X ^ 1 holds, then there exists a constant c depending on x(0) (that is, the 
initial value of x) such that cA™ + s(n) solves R. 

Moreover, in both cases, if the exponential-free part of the canonical form of 
(t)" h(n) is 0, then we can further require that deg(s(n), n) < deg(h(n),n) holds. 

Proof: By Lemma 6, we have 

x(n)=X n ^x(0) + J2^j. (2) 

Denote by terms(/i) all the terms of the canonical form of h(n). Assume each 
t € terms(/i) is of form 

c t n«* ft, 

where c t is a constant in Q, q t is a non-negative integer and j3 t is a prod- 
uct of finitely many elements (with possible repetitions) from {ai, . . . , a/j}. 
Define g(n) := jj^pr- Then g(n) is a poly-geometrical expression in n w.r.t. 
{A}teterms(h). {■ Clearly we have 

g(n)= £ i n *(y)" 

tGterms(/i(n)) 

Therefore, we have 

E at& = E E^^- ( 3 ) 

J=0 teterms(ft) J=0 

According to Lemma 4, for each t € terms(/i), we can find a poly-geometrical 
expression 

«t :=(y)"/tW+a t 
in n over Q w.r.t. ^ satisfying 

i- «t = E;=rofi 9t (f) J ; 

2. ft is a polynomial in n over Q of degree qt ( ii (3t ^ A) or a t + 1 (if /3 t = A), 
and at is a constant in Q; note in the later case, ct n qt (4^)" is a summand 
of the constant term of the canonical form of (j) n h(n) is when viewed as 
a polynomial of the n-exponential variables. 

Therefore, using s t (Vi <E terms(/i)), we can simplify the right hand side of 
Equation (2) to 

lx(0) + ]T a t jA"+ £ / t („)#. (4) 

\ t£ terms f/i) / i£terms(/i) 



Assume for each t € terms(/i), we have f3 t — a^' 1 a^' 2 ■ ■ ■ a^ ,h . 
Define 

f3 t (n) := «) e ^K) et ' 2 ••• 
c := x(0) + ^ a t and s(n) := ^ f t (n)f3 t (n). 

i£terms(/i) ieterms(h) 

It is easy to deduce deg(s(n), a") = max teterms(/l ) (deg(/? t (n), a") < deg(/i(n),a 
Finally, one can easily verify that c and s(n) satisfy the requirements in the con- 
clusion. □ 

Remark 2 In Lemma 7, if X is weakly multiplicatively independent w.r.t. ct\, . . . ,a 
then we know that the exponential-free part of the canonical form of (j) n h(n) 
is 0, without computing the canonical form explicitly. 

2.3 Degree preliminaries 

In this subsection, we review some notions and results on the degree of algebraic 
varieties. Up to our knowledge, Proposition 1 is a new result which provides 
a degree estimate for an ideal of a special shape and which can be applied to 
degree estimate of loop invariant ideals. Throughout this subsection, let K be 
an algebraically closed field. Let F be set of polynomials of K[xi,x 2 , ■ . ■ , x s ]. We 
denote by Vks(F) (or simply by V(F) when no confusion is possible) the zero 
set of the ideal generated by F C K[xi, x 2 , . ■ . , x s ] in K s . 

Definition 5 Let FcK s be an r-dimensional equidimensional algebraic vari- 
ety. The number of points of intersection ofV with an (n—r)- dimensional generic 
linear subspace L C K s is called the degree of V [4], denoted by deg(V). The 
degree of a non- equidimensional variety is defined to be the sum of the degrees 
of its equidimensional components. The degree of an ideal I C K[xi,X2, • ■ ■ ,x s ] 
is defined to be the degree of the variety of I inK s . 

We first review a few well-known lemmas. Note that, for a zero-dimensional 
algebraic variety, the degree is just the number of points in that variety. 

Lemma 8 Let V C K s be an r-dimensional equidimensional algebraic variety of 
degree S. Let L be an (n — r)- dimensional linear subspace. Then, the intersection 
of L and V is either of positive dimensional or consists of no more than S points. 

Lemma 9 Let V C K s be a algebraic variety. Let L be a linear map from K s to 
K k . Then we have deg(L(V)) < deg(V). 

Lemma 10 ([8]) Let I C Q[xi,X2, . . . ,x s ] be a radical ideal of degree S. Then 
there exist finitely many polynomials in Q[xi, X2, . . . , x s ] generating I and such 
that each of this polynomial has total degree less than or equal to S. 



Lemma 11 Let V := W fl| =1 Vi with dim(VK) = r. Then we have 

deg(V) < deg(W) max({deg(K) | i = 1 • • • e}) r . 

Proposition 1 Let X = xi,X2, ■ ■ ■ ,x s and Y — y\, y 2 , ■ ■ ■ , yt be pairwise dif- 
ferent s + t variables. Let M be an ideal in Q[Y] of degree du and dimension 
r. Let fi,f 2 ,--- ! fs be s polynomials in Q[Y], with maximum total degree df. 
Denote by I the ideal (xi — fi,x 2 — f 2 , ■ ■ . ,x s — f s ). Then the ideal J := L + M 
has degree upper bounded by du df . 

Proof: We assume first that M is equidimensional. Let L := l\, l 2 , ■ ■ ■ , l r 
be r linear forms in X, Y such that the intersection of the corresponding r 
hyperplanes and V(J) consists of finitely many points, i.e. Hl := J + (L) is zero 
dimensional. By virtues of Lemma 8, the degree of J equals the maximal degree 
of Hl among all possible choices of linear forms l\, l 2 , ■ ■ ■ , l r satisfying the above 
conditions. 

Let L* := l\, l 2 , . . . , I*, where each I* (j = 1 • • • r) is the polynomial obtained 
by substituting Xi with /j, for i = 1 • • ■ s, in the polynomials lj. Consider the 
ideal L* + M in Q[F]. It is easy to show that the canonical projection map 
Tly onto the space of Y coordinates is a one-one-map between Vc* [M + L*) and 
n Y {V C t+s(H L )). Therefore, V C t(M + L*) is zero dimensional and deg(M + i*) = 
dcg(i?i). Hence, viewing Vc*(M + L*) as 

r 

V C t(M)f) V ct (l*) 

and thanks to Lemma 11, we have deg(Vc*(M + L*)) < dud r ^. Therefore, we 
deduce that deg(J) = maxL deg(M + L*) < du d r j holds, by Lemma 8. 

Assume now that Vc* (M) is not necessarily equidimensional. Let V\ , V 2 , ■ ■ ■ , 14 
be an irredundant equidimensional decomposition of Vc*(M), with correspond- 
ing radical ideals P\, P 2 , . . . , Pk- Then, applying the result proved in the first 
part of the proof to each L + Pi (i = 1 • • • k), we deduce 

deg(J) = Eti deg(7 + P i ) 
<£■=! degCP.K/ 
<Eti deg(P)^. 
= d M d r f , 

where is the dimension of Pi in Q[Y]. This completes the proof. □ 

Remark 3 For J in Proposition 1, a less tight degree bound 

d M d r + s 

can easily be deduced from a generalized form of Bezout's bound, since deg(Vc*+= (M)) 
has degree dM and if of dimension r + s in C t+S . 



Example 4 Consider M := (n 2 — m 3 ), gi := x—n 2 —n — m,g2 := y — n 3 — 3n+l, 
and the ideal J := M + ((71,52)- The ideal M has degree 3, and is of dimension 
1 in Q[n, to]. The degree of J is 9, which can be obtained by computing the 
dimension of 

Q(a,b 7 c,d,e)[x,y 7 m 7 n]/(J + (ax + by + cn + dm + e)), 

where a, 6, c, d, e are indeterminates. The degree bound estimated by Proposition 1 
is 3 x 3 ; which agrees with the true degree. 

3 Invariant ideal of P-solvable recurrences 

In this section, we focus on loops with no branches, where the study of loop 
invariants of such loops reduces to the study of algebraic relations among the 
recurrence variables. In particular, we are interested in those whose assignments 
induce a called P-solvable recurrence. We will first formalize the notion of P- 
solvable recurrence. Then in the rest of this section, we will investigate the shape 
of the closed form solutions of a P-solvable recurrence equation, for studying the 
degree and the dimension of invariant ideal. We will provide degree estimates 
for the the invariant ideal, which is useful for all invariant generation methods 
which need a degree bound, like the proposed polynomial interpolation based 
method and those in [14,15,5]. Last but not least, we will investigate the di- 
mension of the invariant ideal. So that we can get a sufficient for non-trivial 
polynomial invariants of a given P-solvable recurrence to exist. Note that in our 
invariant generation method, we do not need (thus never compute) the closed 
form solutions explicitly. 

A "solvable" recurrence relation is, literally, a recurrence relation which can 
be solved by a closed formula depending only on the index number. The P- 
solvable recurrence relations have poly-geometrical expressions (Definition 1) as 
closed form solutions, which is equivalent to the notion of solvable mapping in [18] 
or solvable loop in [12] in the respective contexts. 

Definition 6 (Univariate P-solvable recurrence) Given a recurrence R : 
x(n+ 1) = Xx(n) + f(n) in K, if f(n) is a poly-geometrical expression in n over 
K, then R is called univariate P-solvable recurrence. 

A multivariate recurrence is called P-solvable recurrence, if the recurrence 
variables can essentially (may need a linear coordinate change) be solved out 
one by one from P-solvable univariate recurrences We can define multivariate 
P-solvable recurrence as follows. 

Definition 7 (P-solvable recurrence) Let m, . . . , nk be positive integers and 
define s := m H hrifc. Let M be a square matrix over Q of order s. We assume 



that M is block- diagonal with the following shape: 



M := 



0n2Xni M n2X „ 2 • 0„ 2Xrlfc 



\0 







n k xn 2 



. M 



n k xn k 



Consider an s-variable recurrence relation R in the variables x ll x 2l 
with the following form: 



,x a and 



fx 1 (n + l)\ 
x 2 (n + 1) 
x 3 (n + 1) 

\x s (n + l)J 



= M x 



/ xi(n)\ 
x 2 (n) 
x 3 (n) 

\x s (n)J 



( flnixl \ 
^2n 2 xl 
f3n 3 X 1 



\ffcn fe xl / 



where fi is a vector of length n\ with coordinates in Q and where i{ is a tuple 

of length ni with coordinates in the polynomial ring Q[x\, . . . ,x ni ^ i-n.-i]; f or 

i = 2,...,k. Then, the recurrence relation R is called P-solvable over Q and the 
matrix M is called the coefficient matrix of R. 

It is known that the solutions to P-solvable recurrences are poly-geometrical 
expressions in n w.r.t. the eigenvalues of the matrix M, see for example [18] . How- 
ever, we need to estimate the "shape" , e.g. the degree of those poly-geometrical 
expression solutions, with the final goal of estimating the "shape" (e.g. degree, 
height, dimension) the invariant ideal. In this paper, we focus on degree and 
dimension estimates. 

We first generalize the result of Lemma 7 to the multi- variable case. 



Proposition 2 Let a\,..., a m e Q \ {1}. Let A G 

matrix in the following Jordan form 

fx ••• 0\ 
1 A ••• 
1 A 



, , — sxs , 

MeQ be a 



•• •• 


\0 



'•. 

A 
1 XJ 



Consider an s-variable recurrence R defined as follows: 

X{n + l) sxl = M sxs X(n) sxl +F(n) sxi, where 
(a) X := x\,X2, ■ ■ ■ ,x s are the recurrence variables; 



(6) F := (/i, f 2 , ...,f s ) is a list of poly-geometrical expression in n w.r.t. a\, . . . , a. 
with maximal total degree d. 

Then we have: 

1. if A = 0, then (/i, fx + f 2 , . . . , fi + f 2 -\ h f s ) solves R. 

2. if A = 1, £/ien i/iere exist s poly- geometric expressions {g\, g 2 , ■ ■ ■ , g s ) in 
a\, . . . , a m such that for each i G 1 • • • s, gi is a poly-geometrical expression 
in n w.r.t. a\, . . . , a m with total degree less or equal than d + i. 

3. if ' X £ {0, 1}, then there exists a solution of R, say (2/1,2/2, • • • , y s ), such that 
for each i = 1, . . . , s we have 

y, := qA™ +g u where (5) 

for each i £ 1 • • • s: (a) c» is a constant depending only on the initial value 
of the recurrence; and (6) gi is like in the case of A = 1. Moreover, assume 
further more that the following conditions hold: 

(i) A is weakly multiplicatively independent w.r.t. a\, . . . , a m ; 
(ii) deg(/j, n) — holds for all j G {1, 2, . . . , s}. 
Then, for all i = 1, . . . , s, we can further choose gi such that deg((/i, n) = 
holds and the total degree of gi is less or equal than max(d, 1). 

Proof: We observe that the recurrence variables of R can be solved one after 
the other, from x\ to x s . When A = 0, the conclusion is easy to verify. The case 
A 7^ is easy to prove by induction on s with Lemma 7. □ 

Proposition 3 Let Ai , . . . , A s , oei , . . . , ce m G Q \{1}. Let M E Q be a matrix 
in the following Jordan form 

/ Ai ••• \ 
e 2 ,i A 2 ••• 
e 3j2 A 3 

'■• '•• '•• '•• 

••• A s _i 
V ••• e s , s _i X S J 

where for i = 2, . . . , s, ti.i-i is either or 1. Consider an s-variable recurrence 
R defined as follows: 

X(n + l) sxl = M sxs X(n) sxl + F(n) sxl , 

where 

1. X := Xi,x 2 , ■ ■ ■ ,x s are the recurrence variables; 

2. F :— (/1, f 2 , . . . , f s ) is a list of poly-geometrical expression in n w.r.t. ai, . . . , a, 
with maximal total degree d. 



Then there exists a solution of R, say (t/i, y 2 , ■ ■ ■ , y s ), such that for each i = 
1, . . . , s we have 

Vi := c t A" + gi, (6) 

where 

(a) Ci is a constant depending only on the initial value of the recurrence and 

(b) gi is a poly- geometrical expression in n w.r.t. 

Ai, . . . , Aj_i, ai, . . . , a m with total degree less or equal than d + i. 

Assume further more that the following conditions hold: 

(i) Ai, A2, . . . , A s is weakly multiplicatively independent; 
(ii) deg(/j, n) — holds for all j G {1,2,..., s}. 

Then, for all i = 1, . . . , s, we can further choose yi such that deg(gi, n) = holds 
and the total degree of gi is less or equal than max(d, 1). 

Proof: We observe that the recurrence variables of R can be solved one 
after the other, from x\ to x s . We proceed by induction on s. The case s = 1 
follows directly from Lemma 7. Assume from now on that s > 1 holds and that 
we have found solutions (y\, 1/2, ■ ■ ■ , y s -i) for the first s — 1 variables satisfying 
the requirements, that is, Relation (6) with (a) and (b). We define 

f(n) = / s (n)-e s , s _i?/ s _i(n + l). (7) 

Note that f(n) is a poly-geometrical expression in n w.r.t. Ai, . . . , A s _i, ai, . . . , a v 
with total degree less than or equal to d+s—1. Moreover, for v e {n, A", . . . , A"_ 1; 
a", . . . , a 1 ^} we have 

deg(/(n),u) < max(deg(/ s (n),w),deg(?/ s _i(n),w)) . (8) 

It remains to solve x s from 

x s (n + l) = X s x s (n) + f(n) (9) 

in order to solve all the variables x\, . . . , x s . Again, by Lemma 7, there exists a 
poly-geometrical expression 

y s :=c s A™ + g s (n), 
where g s (n) is poly-geometrical expression in n w.r.t. 

Ai, . . . , A s _ 1 , ot\, . . . , a m , of total degree upper bounded by d + s. This completes 
the proof of the properties (a) and (b) for y s . 

Now we assume that (i), (ii) hold and we prove the second half of the conclu- 
sion. Observe that we have deg(g s (n), n) = dcg(/(n), n), which is 0, according to 
Relation (8) and the fact that we can choose y s ~i such that deg(y s _i(n), n) = 
holds. Next, we observe that for each 

v e {n, A?,...,A™_ 1; a?,..., <C}, 



we have deg(g s (n),v) = dcg(/(n), v), which is less or equal to deg(t/ s _i(n), v) 
by Relation (8). Therefore, the total degree of g s is less or equal than the total 
degree of y s -i, which is less or equal than max(d, 1) by our induction hypothesis. 
This completes the proof. □ 

Theorem 1 Let R be a P 'solvable recurrence relation. Using the same notations 
M, k, s, F, n\, U2, ■ ■ ■ , nk as in Definition 7. Assume M is in a Jordan form. As- 
sume the eigenvalues Ai, . . . , X s of M (counted with multiplicities) are different 
from 0, 1, with Aj being the i-th diagonal element of M. Assume for each block j 
the total degree of any polynomial in fj (for i = 2 • • • k) is upper bounded by dj . 
For each i, we denote by b(i) the block number of the index i, that is, 

6(i)-l b(i) 

11 J ■ '' ' Y,"r ( 10 ) 

j=i i=i 

Let D\ := ri\ and for allj £ {2, . . . , k} let Dj := dj Dj_\ +rij. Then, there exists 
a solution (t/i, yi, . . . , y s ) for R of the following form: 

yi := CiX? + g h (11) 

for all i £ 1 • • • s, where 

(a) Ci is a constant depending only on the initial value of the recurrence; 

(b) gi is a poly- geometrical expression in n w.r.t. Ai,...,Aj_i, and with total 
degree less or equal than D b ^ . 

Moreover, if {Ai, . . . , A s } is weakly multiplicatively independent, then, for all 
i = l,...,k, we can further choose yi such that deg(gi,n) = holds and the 
total degree of gi is less or equal than Il2<t<b(i) max (<^t> !)• 

Proof: We proceed by induction on the number of blocks, that is, k. The case 
k = 1 follows immediately from Proposition 3. Assume from now on that the 
conclusion holds for a value k = £, with £ > 1 and let us prove that it also holds 
for k = I + 1. We apply the induction hypothesis to solve the first I blocks of 
variables, and suppose that y£ is a solution satisfying the properties in the con- 
clusion. For solving the variables in the (£ + l)-th block, we substitute ye to fi + i 
and obtain a tuple of poly-geometrical expressions in n w.r.t the eigenvalues of 
the first t blocks and with total degree bounded by deDg. Therefore, applying 
again Proposition 3, we can find solutions for the variables in the (£ + l)-th block 
satisfying the properties required in the conclusion. This completes the proof. □ 

Note that the degree estimate in Theorem 1 depends on how the block struc- 
ture of the recurrence is exploited, for example, a 2 x 2 diagonal matrix can 
be viewed as a matrix with a single block or a matrix with two lxl diagonal 
blocks. 

In practice, one might want to decouple the recurrence first, and then study 
the recurrence variable one by one (after a linear coordinate change) to get 



better degree estimates for the poly-geometrical expression solutions, regarded 
as polynomials of n-exponcntial terms as the eigenvalues of the coefficient matrix. 
We will just use a simple example to illustrate this idea. 

Example 5 Consider the recurrence: 

x(n + l)\ /2 0\ fx(n)\ I \ 
y(n + l)\ := 3 x y(n) + x{nf 
z(n + l)J \0 3j \z(n) J \x(n) 3 J 

Viewing the recurrence as two blocks (x) and (y,z), the degree estimate ac- 
cording to Theorem 1 would be bounded by 5 (3 X 1 + 2). 

If we decouple the (y, z) block to the following two recurrences 

y(n + 1) = 3 y(n) + x(n) 2 and z(n + 1) = 3 z(n) + x(nf, 

the we can easily deduce that the degree of the poly-geometrical expression for 
y and z are upper bounded by 2 and 3 respectively, again according to Theorem 1. 

It is easy to generalize the previous results to the case of a matrix M which 
is not in Jordan form. Let Q be a non-singular matrix such that J :— QMQ~ X 
is a Jordan form of M. Let the original recurrence R be 

X(n+1) = MX(n)+F. 

Consider the following recurrence Rq 

Y(n + 1) = JY(n) + QF. 

It is easy to check that if 

(yi(n),y 2 (n),...,y s (n)) 

solves Rq, then 

Q^ 1 {yi(n),y 2 {n), . . . ,y s (n)) 

solves R. Note that an invertible matrix over Q maps a tuple of poly-geometrical 
expressions to another tuple of poly-geometrical expressions; moreover it pre- 
serves the highest degree among the expressions in the tuple. 

We turn now our attention to the question of estimating the degree of the 
invariant ideal of a P-solvable recurrence relation. 

Proposition 4 Let R be an s-variable P-solvable recurrence relation, with re- 
currence variables (xi,x 2 , ■ ■ . , x s ). Let I C Q[xi, x 2 , ■ ■ ■ , x s ] be the invariant 
ideal of R. Denote byl e the extension ofT inQ[xi,x 2 , . . . ,x s ]. Let A — ai,a 2 , . . . 
be the eigenvalues (counted with multiplicities) of the coefficient matrix of R. Let 
M be the multiplicative relation ideal of A associated with variables y\ , . . . , y s . 
Then, there exists a sequence of s poly-geometrical expressions in n w.r.t. eti, a 2 , . . 
say 

/i(n, ajj), ... , f s (n, a?, . . . , aJJ), 



which solves R. Moreover, we have 



X e = (S + M) n Q[x u x 2 ,...,x s ], 

where S is the ideal generated by (x\ — fi(n, yi, . . . , y s ), . . . , x s — f s (n, y\, . . . , y s ) 
in Q[xi,x 2 , ■ ■ ■ ,x s ,n,yi, . . . ,y s ]. 

Proof: The existence of fi, f 2 , ■ ■ ■ , f s follows by Theorem 1 and the fact 
that linear combination of poly-geometrical expressions w.r.t. n are still poly- 
geometrical expressions. The conclusion follows from Lemma 5. □ 

The following lemma is not hard to prove and one can find a proof in [11]. 

Lemma 12 Let R be a P-solvable recurrence relation defining s sequences in 
Q s , with recurrence variables (x±, x 2 , ■ ■ ■ , x s ). Let X be the invariant ideal of R 
in Q[x±, X2, ■ ■ ■ , x s ]; and let I be the invariant ideal of R in Q[x\, x 2 , ■ ■ ■ , x s ] . 
Then X equals to X e , the extension ofX in Q[xi,X2, . . . ,x s ]. 

With Proposition 4 and Proposition 1, we are able to estimate the degree 
of polynomials in a generating system of the invariant ideals. Now we are able 
to estimate the total degree of closed form solutions of a P-solvable recurrence 
without solving the recurrence explicitly. 

Theorem 2 Let R be a P-solvable recurrence relation defining s sequences in 
Q s , with recurrence variables (x\, x 2 , . . . , x s ). Let X C Q[x±, x 2 , . . . , x s ] be the 
invariant ideal of R. Let A = a±,a 2 , . . . ,a s be the eigenvalues (counted with 
multiplicities) of the coefficient matrix of R. Let M be the multiplicative relation 
ideal of A associated with variables y\,...,yk- Let r be the dimension of M. . 
Let fi(n, a™, . . . , a£), . . . , f s (n, a™, . . . , a£) be a sequence of s poly-geometrical 
expressions in n w.r.t. at\, a 2 , . . . , a s that solves R. Suppose R has a k block 
configuration as (m, 1), (n 2 , d 2 ), . . . , (n/j, dk)- Let D\ := m; and for all j e 
{2, . . . , k}, let Dj := dj Dj_\ + nj . Then we have 

deg(X) < dcg(M)Dl +1 . 

Moreover, if the degrees of n in fi (i = 1- ■ ■ s) are 7 then we have 

deg(2)<deg(M)D£. 

Proof: 

Denoting by LI the standard projection from Q to Q : 

(x 1 ,x 2 ,...,x s ,n,y 1 ,...,y s ) i-> (x 1 ,x 2 , . . . ,x s ), 

we deduce by Proposition 4 that 

V(X) = LI(V(S + M)), (12) 

where S is the ideal generated by {x\ — f\{n,y\, . . . ,y s ), . . . , x s -f s (n, yi, . . . ,y s ) 
in Q[xi,x 2 , ■■■ ,x s ,n,yi, ... ,y s ]. 



Thus, by Lemma 9, we have 



deg(T) < deg(S + M). 

It follows from Proposition 1 that 

deg{S + M) <deg{M)D r k +1 , 

since the total degree of /» of R is bounded by according to Theorem 1 and 
the dimension of M is r + 1 is in Q[n, yi, . . . , y s ]. 

With similar arguments, the second part of the conclusion follows from the 
fact that S + M can be viewed as an ideal in in Q[xi, X2, ■ ■ ■ , x s , n, yi, . . . , y s ], 
where M has dimension r. □ 

Indeed, the degree bound in Theorem 2 is "sharp" in the sense that it is 
reached by many of the examples we have considered. Let show two of such 
examples below. 

Example 6 (Example 1 Cont.) The corresponding recurrence only 1 block. 
Denote by A := -v f +1 , v ^ +1 . One can easily check that A is weaklymultiplicatively 
independent. Note the multiplicative relation ideal of A associated with variables 
u,v is generated by u 2 v 2 — 1 and thus has degree 4 and dimension 1 in Q[u,v]. 
Therefore, by Theorem 2, the degree of invariant ideal bounded by Ax l 1 . This 
implies that the degree bound given by Theorem 2 is sharp. 

In the rest of this section, we are going to investigate the dimension of the 
invariant ideal of P-solvable recurrences. This can help to answer the following 
natural question: whether or not the invariant ideal of a P-solvable recurrence 
over Q is the trivial ideal of Q[xi, . . . , x s ]? Note that it is obvious that the 
invariant ideal is not the whole polynomial ring. 

Theorem 3 Using the same notations as in Definition 7. Let Ai, A2, . . . , A s be 
the eigenvalues of M counted with multiplicities. Let M. be the multiplicative re- 
lation ideal of X\, A2, . . . , A s . Let r be the dimension of M.. Let I be the invariant 
ideal of R. Then I is of dimension at most r + 1. Moreover, for generic initial 
values, 

1. the dimension of I is at least r; 

2. if is not an eigenvalue of M and Ai,A2,...,A s is weakly multiplicatively 
independent, then I has dimension r. 

Proof: Assume without loss of genericity that M is in Jordan form. By The- 
orem 1, we deduce that R has a solution (/1, fi, . . . , f s ) as follows 

(ci A" + fci(n), C2 A 2 l + h 2 (n),...,c s A ? s l + h s (n)) , 

where for each i G 1 • • • s, Cj is a constant in Q depending only on the initial value 
of R, and hi is a poly-geometrical expression in n w.r.t. Ai, . . . , Aj_i. Moreover, 
we have 



1. for generic initial values, none of Ci, c 2 , . . . , c s is 0; 

2. if the eigenvalues of M can be ordered in Ai, A 2 , . . . , A s s.t. Ai 7^ 1 and for 
each i e 2 ■ ■ ■ s, Aj is weakly multiplicatively independent w.r.t. A l7 A 2 , . . . , A,_i, 
then we can require that, for all i € 1 • • • s, we have deg(/,, n) = 0. 

Viewing n, A™ (for i = 1, . . . , s) as indeterminates, let us associate coordinate 
variable u a to n, to A" (for i = 1, . . . , s). Denote by V the variety of I in Q 
(with coordinates x\, x 2 , ■ ■ ■ , x s ). Note that we have 

dim(V) = dim(I). 

Denote by W\ , W 2 respectively the variety of M in Q a (with coordinates 

s-f 1 

til, v>2, ■ ■ ■ , u s ) and in Q (with coordinates u , u\, u 2 , ■ ■ ■ , u s ). Note that we 
have 

dim(Wi) = r and dim(W2) = r + 1. 
Consider first the map F defined below: 

F : Q S+1 h+ Q 5+1 
(ito.Ul, • • • ,U S ) -> (Cl til + /1, • • ■ ,c s u s + f s ). 

By Theorem 2, we have V = Fo(W2). Therefore, we have we have dim(I) = 
dim(y) < dim(T^ 2 ) =r + l. 

Now assume the initial value of R is generic, thus we have a ^ 0, for all 
i € 1 • • • s. Let us consider the map Fi defined below: 

Fi : Q s+1 ^ Q s+1 

(tiO, til, . . . ,U S ) -> (tio, Cl til + /l, . . . ,C S ti s + /a)- 

Let us denote by V2 the variety Fi(W2). By virtue of Theorem 2, we have 
dim(V2) = dim(W2) = r + 1. Denote by 77 the standard projection map that 
forgets the first coordin ate, th at is, uo- We observe that V = 77(T^). Therefore, 
we have dim(V r ) > dim(77(V2)) — 1 = r. 

Now we further assume Ai 7^ 1 and for each i e 2 • • • s, Ai is weakly mul- 
tiplicatively independent w.r.t. Ai, A2, . . . , Aj_i the invariant ideal of R. In this 
case, we have that for all z G 1 • • • s, deg(/j,n) = 0. Let us consider the map F2 
defined below: 

(til, ■ • • , U s ) -> (Cl Ml + /l,C2 ti2 + /2, • ■ • ,C S U s + f s ). 

By Theorem 2, we have V — F 2 (Wi). Therefore, we have dim(I) = dim(y) = 
dim(Wi) = r. This completes the proof. □ 

The following result, which is a direct consequence of Theorem 3, can serve 
as a sufficient condition for the invariant ideal to be non-trivial. This condition 
is often satisfied when there are eigenvalues with multiplicities or when and 1 
are among the eigenvalues. 



Corollary 1 Using the same notations as in Theorem 3. If r + 1 < s holds, 
then I is not the zero ideal in Q[xi,X2, ...,x s ]. 

The following corollary indicates that, the fact that the inductive loop invari- 
ant is trivia could be determined by just investigating the multiplicative relation 
among the eigenvalues of the underlying recurrence. 

Corollary 2 Using the same notations as in Theorem 3, consider the corre- 
sponding loop C with xi(0) :— a\, . . . , x s (0) :— a s , where a\,...,a s are inde- 
terminates. If the eigenvalues of R are multiplicatively independent, then the 
inductive invariant ideal of L is the zero ideal in Q[a±, . . . , a s , x\, X2, ■ ■ ■ , x s ]. 

Proof: Since there is only trivial multiplicative relation, the multiplicative 
relation ideal of the eigenvalues is 0, which is of dimension s. By Theorem 3, the 
invariant of R must be zero ideal in Q(di, . . . , a s )[xi, . . . , x s ], since its dimension 
must be at least s. 

Assume there exists a non-zero invariant polynomial p of C, then p must 
be an invariant polynomial of R since the loop variables a\, . . . , a s are free to 
take any value. This is a contradiction to the fact that the invariant ideal of 
R is trivial. Therefore, the inductive invariant ideal of C is the zero ideal in 
Q[ai, . . . ,a s ,xi,x 2 , . • -,x s ]. □ 

Example 7 Consider the recurrence: 

(x(n + l),y(n + l)) := (3 x(n) + y(n), 2 y(n)) with x(0) = a, y(0) = b. 

On one hand, the two eigenvalues are 2 and 3 which are multiplicatively inde- 
pendent, therefore, by Corollary 2, the invariant ideal of the corresponding loop 
is trivial. 

On the other hand, for loop variables (a, b, x, y), the reachable set of the loop 

is 

D\ := {(a, b, (a + b) 3* — 62 l , b2 l ) \ (a,b) G Q 2 , i is a non-negative integer}. 

Therefore, according to Lemma 5, any polynomial vanishes on all points of £H 
must be 0. 

Note in Theorem 3, if we drop the "generic" assumption on the initial values, 
then the conclusion might not hold. The following example illustrate this for the 
case when all the eigenvalues are different and multiplicatively independent, but 
the invariant ideal is not trivial. 

Example 8 Consider the linear recurrence x(n + 1) = 3 x(n) — y(n), y(n + 1) = 
2j/(n) with (x(0),y(0)) = (a,b). The eigenvalues of the coefficient matrix are 
2,3, which are multiplicatively independent. One can check that, when a = b, the 
invariant ideal is generated by x — y. However, generically, that is when a ^ b 
holds, the invariant ideal is the zero ideal. 



4 Algorithm and experimental results 

In this section, we shall discuss how to compute invariant ideals of P-solvable 
recurrences as well as polynomial loop invariants. Our approach is based on 
polynomial interpolation and consists essentially of three main steps. 

1. Sample a list of points S from the trajectory of the recurrence or loop. 

2. Compute all the polynomials vanishing on S up to a certain degree, which 
can be either a known degree bound or a "guessed" bound. 

3. Check whether or not the interpolated polynomials are invariants of the loop. 

As one can see from our algorithm sketch, we need to check whether or not 
a given condition (say a polynomial equation or a polynomial inequality) is an 
invariant. In general, roughly speaking, when a branch condition contains con- 
straints given by inequalities, the problem of checking whether or not a linear 
equation is a loop invariant is undecidable, see [15] for a more detailed discus- 
sion. Nevertheless, criteria showing that a given condition is indeed an invariant 
are useful in practice. For this reason, we are interesting necessary or sufficient 
conditions for a conjunction of polynomial equations to be an invariant of a loop. 

4.1 Checking invariants 

Proposition 5 states a necessary condition for a set of polynomials to be the 
invariant ideal of a given loop. 

Proposition 5 Given a loop C with only one branch and let A be the assignment 
function. Let I be the inductive invariant ideal of C. Then for any point a G 
V(I), we have A(a) € V(I). 

Proof: Denote by T the inductive trajectory of C Let W be the Zariski clo- 
sure of A(V(I)). Let W\ be the Zariski closure of W\V(I). We proceed by contra- 
diction, thus we assume Wi ^ 0. Then we have V(I) = A' 1 {W-^UA' 1 {V {I)f\W) 
and T C A^iVil) n W) and T % A~ 1 {Wi), contradicting the fact that V(I) 
is the Zariski closure of T ■ □ 

The following Proposition, which follows directly from the definition of an 
inductive invariant, can serve as a sufficient condition for a set of polynomials 
to be inductive invariants of a given loop. 

Proposition 6 Let C be a loop with variables X and m branches (Ci,Ai) i = 
1, . . . , m. Let P C Q[A]. IfV(P) contains the initial values of C, and if for each 
a € V(P) n Z(d), we have Ai(a) € V(P), then all the polynomials in P are 
inductive polynomial equation invariants of C. 

Note Proposition 6 states a sufficient condition for a set of polynomials to be 
invariant, not to generate the invariant ideal. 

We shall use Proposition 6 as a criterion to certify given polynomials are 
indeed inductive invariants. Actually, most loop invariant checking criteria work 



in a similar spirit. The proposed criterion is more general than the various "con- 
secutions" conditions in [21], in the sense that all invariants certifiable by those 
"consecutions" conditions is certifiable by the proposed criterion, but there are 
invariants certifiable by Proposition 6, which can not be certified by any of the 
"consecutions" conditions. 

4.2 Implementation of the method 

We use polynomial interpolation to construct candidate invariants from a given 
template (which is either all possible dense polynomials up to a certain degree 
or a specific form guessed by an oracle). To do so, we need to take sufficiently 
many points from the trajectory of the program execution. This is done by 
emulating the program and recording the relevant values. To apply the criterion 
of Proposition 6, we need to compute the image of a variety under a polynomial 
map. This is where we use state-of-art computer algebra software tools. 

In this section, we describe two algorithms for generating polynomial loop 
invariants that we have implemented. We refer the first one as our direct method. 

Notation 3 Notations in the input of our algorithms: 

(i) M :— mi, m2, ■ ■ ■ , m c is a sequence of monomials in the loop variables X 

(ii) S := Si,S2, . . . , s r is a set of r points on the inductive trajectory of the loop 

(Hi) E is a polynomial system defining the loop initial values 

(iv) B is the transitions (Ci,A\), . . . , (C m ,A m ) of the loop 

The subroutines in Algorithm 1 are explained as follows: BuildLinSys(M, S) 
returns anrxc matrix L, such that Lij is the evaluation of the i-th. monomial in 
M at the j-th point in S. LinSolve(L) returns a matrix TV in row echelon form 
with full row rank, whose rows generate the null space of L in Q c . GenPoly(M, 
v) returns the polynomial J2i=i v i m ii where v = (vi,v 2 , ■ ■ ■ , v c ) is a vector in 

Q c . 

Note that we can find effective tools for all the operations in Algorithm 1, 
for instance, we can find tools in [2] for computing the intersection of two con- 
structible sets, or the image of a constructible set under a polynomial map as 
well as testing the inclusion relation. 

However, there is a notable challenge with our direct algorithm (Algorithm 1): 
the coordinates of the points sampled on the trajectory often grow dramatically 
in size. This has clearly a negative impact on the solving of the linear system L. 
All this leads to a severe memory consumption issue, so we decided to consider an 
algorithm based on modular techniques. We opted for a "small prime" approach, 
see Algorithm 2, as we observed that many invariants of practical program loops 
have often small coefficients. 

Some additional subroutines, used in Algorithm 2, are specified hereafter: 
MaxMachinePrimeQ returns the maximum machine-word prime; PrevPrime(p) 
returns the largest prime less than p; BuildLinSysModp(M, S, p) returns an 
r x c matrix L, such that Lij is the evaluation of the i-th monomial in M at 



the j-th point in S modulo p. LinSolveModp(L, p) returns a matrix N in row 
echelon form with full row rank, whose rows generate the null space of L in Z£. 
RatRecon(N, P) returns a matrix N with rational coefficients, such that for each 
i = 1 . . . k, the «-th matrix in N equal to the image of N modulo the i-th prime 
in P if possible; otherwise returns FAIL. 

Proposition 7 Both Algorithms 1 and 2 terminate for all inputs. Moreover, 
when the output is not FAIL, it is a list of polynomial equation invariants for 
the target loop. 

Proof: The termination is easy to check since all loops iterate on finitely many 
terms and each operation (sub-algorithm) does terminate. When the output is 
not FAIL, that means the output satisfies the sufficient conditions for polynomial 
invariants stated in Proposition 6 and thus the conclusion follows from Proposi- 
tion 6. □ 

Remark 4 We handle "unlucky primes" by checking the dimension of the solu- 
tion space (lines 15 — 18 in Algorithm 2): if the dimension of an image increases, 
then we drop this image; if a new image has a lower dimension, then we drop 
all previous images. Several points of Algorithm 2) returns the same 'FAIL" 
message, for sake a simplicity. However, we could customize the FAIL message 
in each return point, for examples: 

— the FAIL at line 9 implies that either the invariant ideal is the zero ideal or 
the total degree of interpolated polynomials is too low; or the modulus is too 
small; 

— the FAIL at line 23 means that the product of the chosen moduli is still too 
small and more images are needed. 



Algorithm 1: PlainInvInterp(M, S, B, E) 
Input: See Notation 3 for M, S, B, E 

Output: A set of polynomial inductive invariants of the target loop 

1 L := BuildLinSys(M, S); 

2 N := LinSolve(L) ; 

3 F:=0; 

4 foreach row vector v € N do 

5 |_ F := F U {GenPoly(T, v)}; 

6 if Z(E) % V(F) then return FAIL; 

7 foreach (d,Ai) e B do 

8 |_ if Ai(V(F)nZ(d)) <Z V(F) then return FAIL; 

9 return F; 



Note that Algorithm 1 and Algorithm 2 will sometimes return FAIL even 
if the bounds for the polynomial degrees and coefficient sizes are known. When 



these algorithms return a list of non-trivial polynomials, we are not sure whether 
those polynomial can generate the whole loop invariant ideal or not. However, 
in practice, these algorithms often find meaningful results quickly. Indeed, both 
algorithms run in singly exponential time w.r.t. number of variables for a fix 
total degree bound, which is stated formally as below. 

Proposition 8 Algorithm 2 runs in singly exponential time w.r.t. number of 
loop variables. 

Proof: The complexity of Algorithm 2 between Lines 1 and 26 is polynomial 
in the number of monomials in the support. 

The number of those monomials is singly exponential t w.r.t. number of loop 
variables. In addition, applying our criterion to certify the result (Lines 27 to 
30) can be reduced to an ideal membership problem, which is singly exponential 
w.r.t. number of loop variables. □ 

In particular, if the total degree bound supplied is greater of equal than 
the degree of invariant ideal and the sample points are sufficiently many, then 
with a high possibility (depending on the selection of sample points and also on 
the choice of the moduli for Algorithm 2), a list of polynomials generating the 
invariant ideal will be computed by out method. 

4.3 Experimental results 

We have applied Algorithm 2 to the example programs used in the paper [20], 
and we are able to find the loop invariants by trying total degree up to 4 for 
most loops within 60 seconds. See the Table 1 for details. 

In the following table, we supply experimental results for computing absolute 
inductive invariants for some well-known programs from literature as well as 
some homemade examples marked with a star *. The first column labeled by 
vars" is the number of loop variables; the second column labeled by "deg" is the 
total degree tried for the methods which use a degree bound; the third column 
labeled by "PI" is the timing of the our method; the fourth column labeled by 
"AI" is the timing of the method described in [5]; the fifth column labeled by 
"IF" is the timing of the method described in [20] ; the sixth column labeled by 
"ALIGATOR" is the timing of the method described in [12]. The time unit is 
second; the "NA" symbol in a time field means that the related method does 
support the input program; the "FAIL" symbol in a time field means that the 
output is not "correct" . All the tests were done using an Intel Core 2 Quad CPU 
2.40GHz with 8.0GB memory. 

The following example shows how we can use the degree and dimension in- 
formatiomw to assure that we are computing the whole invariant ideal. 

1 For more details, see http://www.csd.uwo.ca/~rong/loop_inv.tgz for the source 
of all the programs. 

2 There might be a bug in the version of Aligator we are using, because the compu- 
tation can not finished in lhr in this test; the timing was reported by Laura Kovacs 
in a demo of Aligator. 



Algorithm 2: ModpInvInterp(M, S, B, E, n) 

Input: See Notation 3 for M, S, B, E; n is the maximal number of modular 
images to use 

Output: A set of polynomial inductive invariants of the target loop 

1 p := MaxMachinePrimeQ; 

2 L := BuildLinSysModp(M, S,p); 

3 JV := LinSolveModp(L, p) ; 

4 d := dim(iV) ; 

5 N:= (N); 

6 P := (p); 

7 i := 1; 

8 while i < n and p > 2 do 



9 


if d = then return FAIL; 


10 


JV := RatRecon(N,P); 


11 


if JV / FAIL then break; 


12 


P 


:= PrevPrime(p); 


13 


L 


:= BuildLinSysModp(M, S, p); 


14 


N := LinSolveModp(L, p) ; 


15 


if d > dim(iV) then 


16 




d := dim(AT); N := (N); 


17 




P := (p); 


18 




t:=l; 


19 


else if d — dim(B) then 


20 




N := Append(N, JV); 


21 




P :=Append(P,p); 


22 




i := i + 1; 



23 if i > n or p = 2 then return FAIL; 

24 F ~ 0; 

25 foreach row vector v G JV do 

26 |_ F := F U {GenPoly(T, v)}; 

27 if Z(£) g V(F) then return FAIL; 

28 foreach (Ci, A 4 ) G J? do 

29 |_ if Ai(V(F)nZ(Ci)) <Z V(F) then return FAIL; 

30 return F; 



Table 1. Experiments on selected programs 



prog. 1 


jfc vars 


deg 


_r 1 


A T 
A 1 


r r 




cohencu 


4 


■5 


U.O 


n no 
U.yo 


n oq 


A 1Q 

U.lo 


LUIlfcillL, LI 


4 


2 


u.uu 


u. / u 






format 


5 


4 


3.74 


0.79 


0.37 


0.1 


prodbin 


5 


3 


1.4 


0.74 


0.36 


0.13 


rk07 


6 


3 


3.1 


2.23 


NA 


0.35 


kov08 


3 


3 


0.2 


0.57 


0.22 


0.01 


sum5 


4 


5 


12 


1.60 


2.25 


0.16 2 


wensley2 


3 


3 


0.4 


0.84 


0.39 


0.21 


int-factor 


6 


3 


60.9 


1.28 


160.7 


0.9 


fib(coupled) 


4 


4 


2.4 


0.71 


NA 


NA 


fib (decoupled) 


6 


4 


4.3 


1.28 


160.7 


FAIL 


non-inv2* 


4 


3 


1.2 


3.83 


NA 


FAIL 


couplcd-5-1* 


4 


4 


1.1 


9.58 


NA 


NA 


coupled-5-2* 


5 


4 


5.38 


15.8 


NA 


NA 


mannadiv 


3 


3 


0.1 


0.83 


NA 


0.04 



Example 9 Consider the following recurrence relation on (x,y,z): 

x(n + l)\ /00 1\ /x(n)\ 
y(n + l)\ = 10-3 y(n) 
z(n+l)J \0 1 3 / 

wit/i initial value (x(0),y(0), z(0)) — (1,2,3). Denote by M the coefficient ma- 
trix. Note that the characteristic polynomial of M has 1 as a triple root and the 
multiplicative relation ideal of the eigenvalues is zero- dimensional. So the invari- 
ant ideal of this recurrence has dimension either or 1 . On the other hand, we 
can show that for all k <G N, we have M k ^ M; so there are infinitely many points 
in the set {(x(k),y(k), z(k)) \ k e N} 7 whenever (x(0),y(0), z(0)) ^ (0, 0, 0). 
With our method, we are able to compute the following invariant polynomials 

x + y + z - 6, y 2 + 4yz + 4z 2 - 6y - 24z + 20, 

which generate a prime ideal of dimension 1 (thus the invariant ideal of this 
recurrence), in less than 0.25s. 

5 Concluding remarks 

In this article, we propose a loop invariant computing method based on polyno- 
mial interpolation. We supply a sharp total degree bound for polynomials gen- 
erating the loop invariant of P-solvable recurrences. We supply also sufficient 
conditions for inductive loop invariant to be trivial or non trivial. 



The current implementation is for dense interpolation. However, we observe 
that for loops with sparse polynomials in the assignments, the computed invari- 
ants are often sparse too. As future work, we will investigate suitable sparse 
interpolation techniques for interpolating polynomial loop invariant. 

Acknowledgement. We thank Laura Kovavs and Enric Rodriguez-Carboncll 
for providing the implementations of their invariant generation methods for our 
test. 
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